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We present preliminary results at /3 = 1.95 (a = 0.077 fm) on the first unquenched Nf = 2 + 
1 + 1 lattice computation of the Bk parameter which controls the neutral kaon oscillations in the 
Standard Model. Using Nf = 2 + 1 + 1 maximally twisted sea quarks and Osterwalder-Seiler 
valence quarks we achieve 0(a) improvement and a continuum-like renormalization pattern for 
the four-fermion operator. Our results are extrapolated/interpolated to the physical light/strange 
quark mass but not yet to the continuum limit. The computation of the relevant renormalization 
constants is performed non perturbatively in the RI'-MOM scheme using dedicated simulations 
with Nf = 4 degenerate sea quark flavours produced by the ETM collaboration. 
We getZ?^ G/ (a = 0.077) = 0.747(18), which when compared to our previous unquenched Nf — 2 
determination and most of the existing results, suggests a rather weak B^ GI dependence on the 
number of dynamical flavours. We are at the moment analysing lattice data at two additional j3 
values which will allow us to perform an extrapolation to the continuum limit. 
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1. Introduction 



The mixing K° — Z° plays an important role in the understanding of the physics of CP- 
violation. In the Standard model weak effective Hamiltonian (with three light flavours) at the 
lowest order this oscillation is described only by the operator 

1 



O 



AS=2 



[sr^l-Y 5 )d][sY^l-r 5 )d] 



(1.1) 



The matrix element of the AS = 2 effective hamiltonian can be factorized into a short distance 
contribution which can be computed perturbatively and a nonperturbative long distance contribu- 
tion containing the strong interaction effects. These nonperturbative effects are described by the 
hadronic matrix elements of the renormalized four-fermion operator 



(T?\M> e ff 2 \K°) 
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(Z°|<9 AS=2 Gu)|^°) (1-2) 



The Bk parameter parametrizes the deviation of the hadronic element from the Vacuum Inser- 
tion Approximation (VIA) 



(K°\d AS = 2 (n)\K°) = (K Q \d AS = 2 \K°) VIA B K (ii) 



-f 2 M 2 B K (LL) 



(1.3) 



Combining the lattice computation of the Bk parameter with the experimental value of Sk one 
can constrain the values of the CKM matrix elements. 



2. Lattice setup 

In the gauge sector we use the Iwasaki action while the dynamical quarks have been regular- 
ized employing the twisted mass formalism at maximal twist [|I|] which provides automatic &(a) 
improvement [Q, |[J. The fermionic action for the light doublet in the sea is given by 

s sea,Mtm = £ ^ + ^ + ^ (2 1} 

x 

where we follow the notation in pp. In the heavy sector the sea quark action becomes 

s ,eaMm = £x h (x) [D W [U] +m , h + i^%Tl + Xh (x) (2.2) 

x 

The &(a) improvement and a continuum-like renormalization pattern for the four-fermion 
operators can be achieved by introducing an Osterwalder-Seiler [||] valence quark action allowing 
for a replica of the down (d,d f ) and the strange (s,s') quarks The action for each OS valence 
flavour %f reads 

s val,OS = £ ^ {u] + ^ + iwrf] ^ (2 3) 

x 

where the Wilson parameters should satisfy the relation —r s = r c j = = r s i. 

In table [j] we give the details of the simulation of the bare Bk parameter. For the inversions 
in the valence sector we used the stochastic method with propagator sources located at random 
timeslices in order to increase the statistical information [|7| ^J. 
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a\i sea 


aM% (MeV) 


L 3 x T 


#confs 




0.0025 


~270 


32 3 x 64 


144 


j3=1.95 


0.0035 


~320 


32 3 x 64 


144 


(a ~ 0.077fm) 


0.0055 


~400 


32 3 x 64 


144 


a\i a = 0.135 a/J-s = 0.17 


0.0075 


~460 


32 3 x 64 


80 


a\i wl = {aAt sea ,0.0141,0.0180,0.0219} 


0.0085 


~490 


24 3 x 48 


244 



Table 1: Simulation details for the bare Bk parameter. The bare valence quark mass equal to 
0.0141,0.0180,0.0219 have been chosen to allow for a smooth interpolation around the physical 
strange quark mass. 



ensemble 


a\l sea 


am PCAC 


aM sea 


Qsea 


a\l val 


am PCAC 


#confs 


lm 


0.0085 


-0.04125(13) 


0.03308(10) 


-1.3109(08) 


set 1 


-0.02116(2) 


304 


IP 


0.0085 


+0.04249(13) 


0.03286(09) 


+1.3091(08) 


set 1 


+0.01947(19) 


304 


7m 


0.0085 


-0.03530(13) 


0.02851(10) 


-1.2681(10) 


set 1 


-0.0216(2) 


400 


7p 


0.0085 


+0.03608(11) 


0.02854(08) 


+ 1.2683(09) 


set 1 


+0.01947(19) 


400 


8m 


0.0020 


-0.03627(11) 


0.02804(08) 


-1.4994(02) 


set 1 


-0.0216(2) 


400 


8p 


0.0020 


+0.03624(13) 


0.02743(10) 


+1.4978(03) 


set 1 


+0.01947(19) 


400 


3m 


0.0180 


-0.0160(2) 


0.02191(09) 


-0.6068(59) 


set 2 


-0.0160(2) 


352 


3p 


0.0180 


+0.0163(2) 


0.02183(09) 


+0.6015(57) 


set 2 


+0.0162(2) 


352 


2m 


0.0085 


-0.02091(16) 


0.01815(11) 


-1.0834(32) 


set 1 


-0.0213(2) 


352 


2p 


0.0085 


+0.0191(2) 


0.01692(13) 


+1.0445(45) 


set 1 


+0.01909(18) 


352 


4m 


0.0085 


-0.01459(13) 


0.01404(08) 


-0.9206(43) 


set 2 


-0.01459(13) 


224 


4p 


0.0085 


+0.0151(2) 


0.01420(12) 


+0.9289(64) 


set 2 


+0.0151(2) 


224 



Table 2: Details of the analysed ensembles for the RCs at j8 = 1.95; set 1 = {0.0085, 0.0150, 
0.0203, 0.0252, 0.0298} and set 2 = {0.0060, 0.0085, 0.0120, 0.0150, 0.0180, 0.0203, 0.0252, 
0.0298}. Here M sea = yj (Z A m s ™ AC ) 2 + (n sea ) 2 , with a self-consistent Z A -value. 



The computation of the renormalization constants (RCs) for the relevant two- and four-fermion 
operators has been performed adopting the RF-MOM scheme These RCs are computed by ex- 
trapolating to the chiral limit the RCs estimators measured at several quark mass values. We have 
performed dedicated runs with Nf = 4 degenerate sea quarks in order to be able to estimate the 
chiral limit of the RCs. In these Nf = 4 simulations working at maximal twist would imply a con- 
siderable fine tuning effort due to the difficulties in determining ampcAC near atnpcAC = 0- Instead, 
working out of maximal twist the stability of the simulations increases and the &(a) improvement 
of the RC estimators is achieved by averaging simulations with an equal value of the polar mass 

Za m sea 1773 

M sea but opposite value of m p e ^ AC and , where tan0 = — ^ AC [|10Q. We label these ensem- 
bles as Ep/m where E=l,2... and p/m refers to the sign(0 VfM ). In table || we report the parameters 
of the ensembles analysed to compute the RCs. 
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Figure 1 : (a) Effective mass of the pseudoscalar mesons made out of two mass-degenerate quarks 
regularized with opposite Wilson parameters (see [ |l2| ] for definitions) (b) Quality of the ratio 
plateaus at fixed value of = 0.0180 



3. The K-meson bag parameter 

The only four-fermion operator relevant for Bk in the Standard Model is 1 (see [0, 1 1 ] for 



notation) 

Qi" = °vv+aa = [(sY^)(ffd r ) + {s yii y 5 d){s'y^y 5 d')\ + [d +> d'} (3.1) 

It is multiplicatively renormalizable if the Wilson parameters appearing in the action of the OS 
fermions satisfy the relation — r s = rj = = ry [Hj. 

In order to estimate the ^-parameter we compute a three point correlation function with 
an insertion of free to move in lattice time t, and two 'X-meson walls" with pseudoescalar 
quantum numbers and with a fixed time separation between them \tt — tg\ = T/2. In this way, the 
lattice estimator for the bare Bk parameter can be obtained from the ratio Ri (?) for times ?/. <C t <C ?r 

cf(t-t L )cP (t-t R ) fe«<* 3 

where the two- and three-point correlation functions are defined as in [ 12]. In figures [Ta| and lb the 
quality of the plateaus of the effective mass and the ratio R\ are shown for two values of the light 
quark mass and one typical value of \ih = 0.0180. 

Chiral extrapolations in the light quark mass are performed using SU(2) Partially Quenched 
Chiral Perturbation Theory at NLO JT3| , p4| j. We carry out three chiral fits, one for each simulated 
heavy mass, according to the fit ansatz 



B K (M lh M hh ) =B K {M hh ) 



(3.3) 



as we illustrate in figure ga[ The decay constant fo has been obtained by analysing our data for the 
mass and decay constant of the // pseudoscalar meson following the procedure of [ft]. Our result is 



1 Note that it contains the same physical information as |TTT| 
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fo = 120.99(09) MeV. In order to estimate the systematic errors affecting this extrapolation in the 
w/J-quark mass, we also tried a first order polynomial fit. 

Since we have obtained our lattice data for masses around the physical strange quark mass, we 
determine the bare Bk parameter through an interpolation in (Mhh) 2 to the physical mass M% h = 
2M K — M 2 as is shown in figure 2b. 

Finally, the renormalized Bk parameter is given by the relation [|12|], 

Z+ 

B K = —^rB K (3.4) 

where Zj'j, and Zy are the RCs of the operators Q^, the axial and the vector currenty. The 
relevant RCs both for the two- and four-fermion operators have been computed non-perturbatively 
in the RF-MOM scheme following Our strategy to remove the inescapable G(a) discretization 
effects consists in performing a 9 sea -average of the RCs estimators as described in Q10|]. 

Using the notation of flTij], we compute the elements of the dynamical matrix D i; for each 
value of a 2 p 2 = £ v sin 2 (ap v ), \i val and Ep/m (and hence of M sea ). As described in jl2|], the 
valence chiral limit extrapolation is safely determined with the help of the following fit ansatz 

D u (p 2 , Ep/m;^ 1 ) = D, 7 (^;Ep/m) + A ;7 (/;Ep/m)M™' + ^^'^ (3.5) 

\ PS ) 

at fixed values of a 2 p 2 and Ep/m (and hence of M sea ). The last term takes into account the con- 
tribution of one Goldstone boson pole (GBP) (see Appendix A of [12] for details). In figure |a| 
we show the smooth dependence of D\\ on {aM^) 2 for a representative value of the momentum 
even before the GBP substraction. Hence the one GBP contribution is strongly suppressed on the 
vertex D\\. Having performed the GBP subtraction and the valence chiral extrapolation, the 6{a) 
discretization effects are removed by averaging D ;y (^ 2 ;Ep/m; i u V0 ' = 0) over 6 sea . By combining 
the valence chiral limit estimator of Dn with the corresponding chiral limit of the quark field RC, 
Z q , we obtain the intermediate quantities Z^(p 2 , pL val = 0,Ep/m). Then, the sea chiral limit is taken 
at each fixed value of (ap) 2 by fitting these intermediate estimators to a polynomial of first order in 
(aM sea ) 2 . We find that the dependence on the sea quark mass is very mild, as it can be seen from 
figure |3^. In this way, we get the chiral limit RC estimators, ZJj (p 2 ). 
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(aft 2 (aM sca ) 2 
(a) (b) 

Figure 3: (a) Goldstone boson pole substraction and valence chiral limit of D\\ for the ensemble 2p, 
ix sea =0.0085 and (a/5) 2 =1.889104 ; (b) Sea chiral limit extrapolation of Z+ for (a/5) 2 = 1.889104. 




(ap) (ap)' 
(a) (b) 



Figure 4: (a) Mixing coefficients in the valence and sea chiral limit ; (b) RI'-MOM computation of 
the multiplicative renormalization factor Zj^ . 



An important issue that arises in our analysis of the RCs is that wrong chirality mixings A,-,- 
can affect the renormalization pattern of at order a 2 or higher [Q]. However, as it is clearly seen 
from figure ^a[ these mixing coefficients are negligible within errors in our data. 

Improved chiral limit RC estimators are obtained by analytically subtracting 6(a 2 g 2 ) dis- 
cretization errors directly from the vertex Djj, using the one-loop results in lattice perturbation 
theory calculated in [15]. The effect of this perturbative correction is illustrated in figure ^ where 
we compare the RI'-MOM values for Z^ at the reference scale /i = a~ l in three cases: uncor- 
rected, corrected with a bare lattice coupling go = 6/j8 2 and with a boosted one g 2 h = g^/(P), 
where the average plaquette (P) is computed non perturbatively. 

After bringing Z^ (p 2 ) to a common reference scale p 2 = ji^ = a~ 2 by employing the known 
NLO running formula [jig), we applied two different methods to remove the remaining ff(a 2 p 2 ) 
discretization errors [|TJ. In the Ml-method we perform a linear fit in (ap) 2 of ZJj(/1q) i n the 
interval (ap) 2 = [1.5 : 2.0] (see figure 4b). Alternatively, the M2-method consists in simply taking 
the weighted average of the RC in the momentum window (a p) 2 = [1.8 : 2.0]. 

Combining four- and two-fermion RCs computed as in [p^|], we obtain our final result for the 
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RGI Bk at j8 = 1.95, corresponding to a = 0.077 fm, by employing the estimates of the RCs from 



the procedure Ml and eq.(p.3|) for the extrapolation of our Bk data to the u/d-quark physical point 

B K GI (a = 0.077) = 0.747(05)(17)[18] 
where the first error is statistical, estimated from a boostrap analysis, and the second results from 



summing in quadrature several systematic uncertainties (0.014 from the spread in using eq.(p.3|) 
or a polynomial fit in the extrapolation to the physical point and 0.009 associated to the renormal- 
ization). The error quoted in brackets is the total error obtained by a sum in quadrature of the 
systematic and statistical ones. 
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